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Abstract 

The magnetic field integral equation for axially symmetric cavi¬ 
ties with perfectly conducting piecewise smooth surfaces is discretized 
according to a high-order convergent Fourier-Nystrom scheme. The 
resulting solver is used to accurately determine eigenwavenumbers and 
normalized electric eigenfields in the entire computational domain. 


1 Introduction 

This work is on a numerical solver for the time harmonic Maxwell equations 
in axially symmetric microwave cavities with piecewise smooth and perfectly 
electric conducting (PEC) surfaces. We use the interior magnetic field in¬ 
tegral equation (MFIE) together with a charge integral equation (ChlE) 
and high-order convergent Fourier-Nystrom discretization to find normal¬ 
ized electric eigenfields to high accuracy. 

The intended primary application of our solver is in computational ac¬ 
celerator technology. Our experience is that our solver more than doubles 
the range of frequencies for which electric and magnetic eigenfields can be 
accurately evaluated, in comparison with finite element method programs 
commonly used. This opens up for improved evaluation of, so called, wake- 
fields. Wakefields affect particle trajectories during accelerator operation 
and Wakefield prediction is therefore of great importance in accelerator de¬ 
sign, see m Chapter 11]. 

Our solver is based on the work [20], which in turn draws on progress 
in A major step forward from (201 is the efficient treatment 

of piecewise smooth surfaces and field singularities at sharp edges. For this, 
we rely on a method called recursively compressed inverse preconditioning 
(RCIP) |2I]. 
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The RCIP method can be seen as an automated tool to enhance the 
performance of panel-based Nystrom discretization schemes for Fredholm 
second kind integral equations in the presence of boundary singularities. 
For the determination of normalized eigenfields in microwave cavities with 
sharp edges, and from a numerical point of view, it is important to resolve 
boundary singularities and their associated non-smooth fields to high preci¬ 
sion. This is so since these fields give non-negligible contributions to electric 
and magnetic energies needed in the normalization. From a more practi¬ 
cal point of view in the accelerator design process, the identification and 
evaluation of field singularities is necessary since strong fields may cause 
field emission and quenching (thermal breakdown) in superconducting cav¬ 
ities j26l Chapter 11 and 12]. 

A common approach to the numerical resolution of fields at sharp edges 
is to exploit a priori knowledge of asymptotic behavior and to include a 
leading order singularity, or multiple non-integer powers, in tailor-made ba¬ 
sis functions. This approach generally reduces the convergence order due to 
a dense spacing of presumptive exponents [5j and is difficult to automate 
and to apply to problems that are not translationally invariant in one direc¬ 
tion [22]. Still, it has been used in numerous papers where the method of 
moments (MoM) is applied to scattering from PEC structures with sharp 
edges |5] and also to find stress fields around cracks, notches, and grain 
boundary junctions in computational mechanics [23]. The RCIP method, 
on the other hand, does not require any known asymptotics and generally 
retains the convergence order of the underlying discretization scheme. 

The construction of our solver covers a wide range of topics and com¬ 
putational techniques that are more or less well known and it would carry 
too far to review them all. Some techniques that are considered particularly 
important are discussed as they appear in the text. For other issues we 
merely give references. 

The outline of the paper is as follows: Section presents the MFIE 
and the ChlE for the problem at hand and an integral representation for 
the electric field in a concise notation. Section Him and review the 
Fourier-Nystrom discretization scheme for smooth surfaces. The emphasis 
is on kernel evaluation and on the conversion of a volume integral, used 
for normalization, into an expression that is more suitable for numerics. 
Section is about sharp edges and how the RCIP method is incorporated 
into the scheme. Here the ChlE plays an important role by simplifying the 
accurate extraction of the surface charge density. Section relates the com¬ 
puted complex valued electric fields to physical time-domain standing wave 
fields. Section contains numerical examples with relevance to accelerator 
technology and Section 10 discusses future work. In order to maintain a 
high narrative pace in the main body of the paper, and also to provide an 
overview and to facilitate coding, all explicit information on the integral 
operators used is gathered in an appendix. 
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Figure 1: An axisymmetric surface F generated by a curve 7. (a) A point r on 
F has outward unit normal v and tangential vectors r and 0. (b) r has radial 
distance p, azimuthal angle 6, and height z. The planar domain A is bounded 
by 7 and the 2 ;-axis. (c) Coordinate axes and vectors in the half-plane 0 = 0. 


2 Problem formulation 

This section introduces the MFIE for the time harmonic Maxwell equa¬ 
tions in a notation that is particularly adapted to electric fields inside ax¬ 
ially symmetric cavities with PEC surfaces. Parts of the material are well 
known [nmansiEb]. The presentation parallels Section II of [2^ , in which 
magnetic fields are of primary interest. 

2.1 Geometry and unit vectors 

Let F be an axially symmetric surface enclosing a three-dimensional domain 
V (a body of revolution) and let 

r = (x, y, z) = (pcosff, psinff, z) (1) 

denote a point in M^. Here p = is the distance from the z-axis 

and 0 is the azimuthal angle. The outward unit normal iz at a point r on F 
is 

u = {vp cos 9, i^p sin 9, Vz) ■ (2) 

We also need the unit vectors 

p = (cos 9, sin 9, 0), 

0 = (—sin0, cos 0,0), ^ ^ 

(3) 

T = 6 X V = {vz cos 0, Uz sin 0, —Vp ), 

^ = ( 0 , 0 , 1 ), 

where 0 and r are tangential unit vectors. See Figure [^a) and[2b). 
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The angle 0 = 0 defines a half-plane in whose intersection with T 
corresponds to a generating curve 7 . Let r = (p, z) be a point in this half¬ 
plane and let A be the planar domain bounded by 7 and the z-axis. The 
outward unit normal on 7 is = (up, v^) and r = —i^p) is a tangent. See 

Figure [^c). The unit vectors in the p- and z-directions are p and z. 

2.2 PDE formulation 

The electric field is scaled with the free space impedance 70 such that E[r) = 
PQ^E'{r), where E'{r) is the unsealed field. With vacuum in V and with 
r perfectly conducting, the electric field E{r) satisfies the system of partial 
differential equations 

V^E{r) + k^E{r) = 0, r£V, (4) 

V-E{r) = 0, r G 1/, (5) 

with boundary condition 

lim 1 /° X E{r) = 0, r° G T. ( 6 ) 

V Br^r° 

We will find nontrivial solutions to these equations in a fast and accurate 
fashion via the MFIE. Note that, since F is perfectly conducting and from 
a mathematical as well as physical point of view, nothing exterior to F can 
affect E{r) in V. Hence the region exterior to F is irrelevant to the problem. 

The values k'^ for which the system Q, ([^, and Q admits nontrivial 
solutions are called eigenvalues. We refer to the corresponding fields E{r) 
as electric eigenfields and to k as eigenwavenumbers. The eigenvalues con¬ 
stitute a real, positive, and countable set, accumulating only at infinity m- 
The eigenvalues have finite multiplicity. Electric eigenfields that correspond 
to distinct eigenvalues are orthogonal with respect to the inner product 

{E,G)= [ E*{r)-Gir)dV, (7) 

Jv 

where E{r) and G{r) are vector fields on V and the asterisk indicates the 
complex conjugate. Subspaces of electric eigenfields that correspond to de¬ 
generate eigenvalues can be given orthogonal bases. 

Electric eigenfields E(r) are normalized so that 

\\Ef= [ E*{r) ■ E{r)dV = 1. (8) 

Jv 

The volume integral in Q is referred to as the normalization integral. In 
Section]^ it is converted into a surface integral that is well suited for numer¬ 
ical evaluation in the framework of the MFIE. 
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2.3 Integral representation of the electric field 

The surface current density Jg and the surface charge density Qs are defined 
as 


JJr°) = I lim i/° X (V X E(r)), r° e T , (9) 

k V'9r-s>r° 

Qs{r°) = lim —u°-E{r), r° G T. (10) 

V 9r—>r° 


When Jg and are known the electric field is given by the integral repre¬ 
sentation 


E{r) 

0 


ik J Jg(r')^>fc(r,r') dr' 


Pg(r')V$fc(r,r')dr' 


r€V, 
r € cV, 
( 11 ) 


where cV is the exterior to 1/ U T. Here, using the time dependence e 
with angular frequency a; > 0, the kernel 




giA:|i—r'l 

47 r|r — r' 


( 12 ) 


is the causal fundamental solution to the Helmholtz equation and 

V^k{r,r') = —^gPdr-r'l), (13) 

47r|r — r'\-^ 

where 

P{\r - r'l) = (1 - \k\r - r'|)e'^l'’-^'l . (14) 

The lower equation in 0 states that Jg and induce an electric null held 
outside r. 

We decompose the electric held in its cylindrical coordinate components 


E{r) = pEp{r) 0Eg{r) + zE^{r) 
and the surface current density in its tangential components 


(15) 


Js{r) = TJr{r) + eje{r). (16) 

From ( |11[ ) the components of the electric held, including the induced external 
electric null held, can be expressed as 

Ep{r) = ikSjJrir) + kSsJe{r) + KiiQs{r ), 

E0{r) = kSgJrir) + ikSioJeir) + ii^i2&(r), (17) 

E^{r) = ikSuJrir) + Ki3Qs{r) , 
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where r G \ F and the double-layer type operators Kq and single-layer 
type operators Sa with various indices a are defined by their actions on a 
layer density g{r\ r G F, as 


and 


Kag{r) 

= J^Ka{r,r')gir')dT', 

(18) 

Ka{r,r') 

= Dair,r')P{\r -r'\), 

(19) 

Sag{r) 

= J^SUr,r')gir')dr', 

(20) 

Sa{r,r') 

= Z„(r,r')e‘^l^-^'l. 

(21) 


The functions Da{r,r') and Za{r,r') can be viewed as static kernels, cor¬ 
responding to wavenumber k = 0. Their explicit expressions are given in 


Appendix A.l 


2.4 Integral equations for Jg and gg 

The interior MFIE for Jg reads 

J,(r) - 21/X ^ (J,(r') X V4>fc(r,r')) dF' = 0, rGF. (22) 

The surface charge density Qg is related to the surface current density Jg by 
the continuity equation 

Gsif") = ‘ , (23) 

where Vs-( ) is the surface divergence. As in |20j . we avoid the differentiation 
inherent in (|23[) by evaluating Qg from the interior charge integral equation 


Qg{r)-2 j V ■V^k{r,r')Qg{r')(lT' = 


-2ik j u ■ Jg{r')^k{r,r')dT', rGF. (24) 

When Jg and Qg are solutions to (22) and (24), then E{r) in ( [IT| ) is a solution 
to the system (Q, 0, and 0. 

The Fredholm second kind integral equation (24) and its exterior coun¬ 


terpart are here denoted the ChIF. Over the last decade, the ChlF has 
become a tool for dealing with certain numerical problems known as “low- 
frequency breakdown” and which occur when the MFIE, or the related for¬ 
mulations EFIE and CFIE, are used for exterior electromagnetic scattering 
at low frequencies [29]. Low-frequency breakdown is caused by decoupling 
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of electric and magnetic fields. It manifests itself when integral equations 
are solved or when fields are reconstructed. More generally, the ChlE has 
been combined with the MFIE EH]) with the EEIE [Hj, and with the 
CFIE [H |28| for PEC surfaces and with the EFIE and the MFIE for pene¬ 
trable objects [28|. The combination of the MFIE and the ChlE for exterior 
problems is denoted the ECCIE in [29|. Low-frequency breakdown is not an 
issue when computing eigenfields since cavity resonances are wave phenom¬ 
ena with a strong coupling between electric and magnetic fields. Our reasons 
for preferring the ChlE over (23) have, as in [20], to do with convergence 
order and achievable accuracy. 


In accordance with m we split ([ 2 ^ into the two coupled scalar equations 


(I + K,) Mr) + iK2Jeir) = 0, 
iKsMr) + (I + M) Je{r) = 0. 


(25) 


and write (24) as 

(/ - 2M)Q,{r) = -2ikS5Jrir) + 2kSeJe{r). 


(26) 


Here I is the identity. The operators Ka, a = 1,2, 3,4, v are of the double¬ 


layer type (18) with (19). The operators ^5 and Se are of the single-layer 


type (20) with (21). See Appendix A.l for their explicit expressions. 


3 Fourier series expansions 


The aim of this paper is to present a high-order convergent and accurate 
discretization scheme to solve the MFIE and to evaluate electric eigenhelds, 
normalized by ([^. We employ a Fourier-Nystrom technique where the first 
step is an azimuthal Fourier transformation of the MFIE system (25) and 
(26) and of the system for the decomposed electric field 

We define the azimuthal Fourier coefficients for 27r-periodic functions 


9nir) 
Gn{r, r') 


1 


1 


/ 71 

-1 


e-'^^g{r)A9, 

-e'), 


(27) 

(28) 


where g{r) represents functions like Qsir), Jrir), Je{r), Ep{r), EqM, and 
Ez{r) and where G{r, r') represents Ka{r, r'), Da{r, r'), Sa{r, r'), Za{r, r'), 
and P{\r — r'\). The subscript n is the azimuthal index, n = 0, ±1, ±2,.... 

We also define the modal integral operators K^n and San in terms of 
their corresponding Fourier coefficients Kan{r,r') and San{r,r') as 

Kangnir) = V^ f Kan{r,r')gn{r')pd'y', (29) 

J7 
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(30) 



r, r')gnir')p d'y’. 


Expansion and integration of (25) over 9' now give the system 




integral equations 


(/ + Kin) Jrnir) + iK 2 nJen{r) = 0 , 
iKsnJrnir) + {I + Kin) Jdn{r) = 0 , 


where r G 7 . Analogously, the modal version of (26) is 


(/ 2Ki/n) £*sn(^) — ‘^^kSllnJrnij') ‘^kSflnJonij') • 


The modal representation of the electric field in (17) is 


— diS'inJTnij’) T kSgnJoni'^) T KiinQsn{,'l") j 
^Bnir) — kSQnJrni^) T ^kSiQnJoni^) ^Ki2nQsni'l") ; 
k^znij") — 'diSiinJrnij') T Ki^nQsni^) ; 


(31) 


(32) 


(33) 


where r ^ 7 . 

In what follows, we sometimes present only modal expressions for opera¬ 
tors and helds. The reasons being the close resemblance between expressions 
in and modal expressions, visible in the examples above, and that we seek 
eigenfields 


En{r) = (pEpnir) + OEonir) + zE^nir)) 


'mO 




for one azimuthal mode at a time. 


(34) 


4 Conversion of the normalization integral 

The normalization (volume) integral in ([^ can be converted into a sum of 
line integrals over 7 that is more suitable for numerical evaluation. The 
line integrals contain Fourier coefficients of an electric scalar potential 'l'(r) 
and a magnetic vector potential A{r). In [20l Appendix A], using results 
from [3], this conversion is done in the context of magnetic eigenfields. Since 
for each eigenwavenumber the electric and magnetic eigenfield energies are 
equal, the expression in [20] applies equally well to electric eigenfields. With 


En (r) as in (|34[) the converted expression reads 


\Er 


= K 


[ {\An\^ -\^n\^)pdA+ [ {A*n-J,n + 'lkAln'i>n)pd^, (35) 
Ja J-y 
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where 


k'^ [ \'^n\‘^pdA = -I- f {n^ - k‘^p‘^)\^n\‘^ dj 
Ja ^J-yP 

+ ^ Re {( 2 t • r[dr^*)n + ^;) {d,+ ^)n] pd^ , ( 36 ) 

k"^ f \An\^pdA = -\ I[{n^ - k"^p^)\An\^) d7 


'1 


P 


1 f u ■ r 


(|Ap„p + - 4nIm{A*„Ae„}) d 7 

^ Jy h- 

~\j ’^''^{\i9-r^)n\^-\id,,+ A)n\^'^pd'y 
+ ^e{{‘^T ■ r{d-rA*)n + A*J ■ {d^+A)n} pdj, (37) 

and where directional derivatives of a fnnction g{r) are abbreviated as 


drgir) = T ■Vg{r), reT, 
d^+g{r°) = lim u° ■ Vg{r ), r° G T . 

VBr^r° 


(38) 


In order to evaluate ( |35| ) from the solution to (31) and (32), the Fourier 
coefficients of 'l'(r) and A{r), and their derivatives with respect to r and 
, need to be related to gsn{r), Jmir), and Jen{i")- Partial information on 
these relations, along with derivations, can be found in |20l Appendix B], 
We now give complete information without derivations. 

Two different decompositions of A„(r) are needed 


An(r) = TArn{r) + 6Aon{r) + uAunir ), 
An{r) = pApn{r) + OAen{r) + zAzn{r). 

Expressions for An(r) are in these bases 

Arn(^) — SinJrni'l"') T '^S2nJ9n{^) j 
AoniX^ — '^^SnJrni^') T S/^nJoniX^ ) 
Aeo{r) = 0, 

Ai/niX^ — S^nJrniX^ A '^S^nJoniX^ ; 
Apn(^) — S\2nJTn{x') T ^^ISnJoni.^') j 
AzniX^ — ^UnJrni'f') ■ 


(39) 


( 40 ) 


The modal operators San, a = 1,..., 6,12,13,14, are defined via (30), (28), 


(21) and Appendix A.l 
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Expressions for the Fourier coefficients of the cylindrical components of 
the derivatives of A(r) with respect to r and are 

~ i-^16n«^Tn(^) “1“ Kx'jnJdnij’^ t 
{d-rA^ — -^18n«^Tn(^) ; 

(^1/+■^)pn(^) — 2 ^2'^Tn(^) “1“-^19n«^Tn(^) “1“ i'^20ri>^l9n(^) ) (^1) 

(^j/+-^)0n(^) — 2^'^21n'^Tn(^) “1“-^22n'^(9n(^) ) 

{di,+A')zn{'l') — “1“ -^23n'^™(^) ■ 


The modal operators Kan, oi = 14, ...,23, are defined via (29), (28), (19) 
and Appendix |A.1[ 

Expressions for the Fourier coefficients of T(r) and its derivatives are 

k 


'^nir) = -pAenir) 
n 


n / 0, 

'^o{r) = S^oQso{r) , 

{d.r^)n{r) = \kArn{r) , 

{d^+^)n{r) = Qsnir) + ikA^nir), 


(42) 


with S^o defined via ([^, ([^, ([^ and Appendix |A.l 


5 Fourier coefficients of static kernels in analytic 
form 

When r and r' are far from each other, all kernels Ka{r,r') and Sa{r,r') 
are smooth functions of 6 — 0' and we evaluate the corresponding Fourier 


coefficients Kan{i’,r') and Sanif',r'), needed in (31), (32), ([3^ and (35), 


from (28) using discrete Fourier transform techniques (FFT). When r and 
r' are close, the kernels vary more rapidly and FFT alone is not efficient. 
Instead we split each Ka{r, r') and ^^(r*, r') into two parts: a smooth part, 
which is transformed via FFT, and a rapidly varying part, which is trans¬ 
formed by convolution of Dcm{r,r') and Zani'f',r') with parts of Pn{r,r'). 
See [ini Section 6] for details on this splitting and [T9l Section 12.1] for a 
definition of when r and r' are considered close. 

The coefficients Dan{r,r') and Zanif^r') can be expressed in terms of 
half-integer degree Legendre functions of the second kind m Equation 
8.713.1] 

cos(nt) dt 


2 J-^ y/8 {x - cos{t)) 


(43) 
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We use these analytic expressions in the convolutions where, in our setting, 


X = 1 + 



(44) 


The functions n„_i(x), with real arguments X > 1 and often with cosh(x) 
substituted for X) may also be called toroidal functions [T6l Section 8.850], 
[211 page 201], ring functions [H Section 8.11], or toroidal harmonics [TB] . 
They are symmetric with respect to n, exhibit logarithmic singularities at 
X = 1, and are relatively cheap to evaluate. 

The particular combinations of toroidal functions 


^n(x) = (x^^n-i (x) - (x)) (45) 


play an important role in our analytic expressions. They multiply functions 
which may exhibit Cauchy-type singularities at x = 1- The 91„(x) are finite 
at X = 1) but have logarithmic singularities in their first (right) derivatives. 

The toroidal functions can be evaluated via a recursion whose forward 
form is 

4-77 _ A 977 _ 

i3n-i(x) = , n = 2, 3,... . (46) 


When 1 < X < 1-0005 we use (46) as it stands. When x > 1.0005, and for 
stability reasons, we use a backward form of (46). A short Matlab code 
which evaluates £}_i(x) and Oi(x) is presented in [20l Appendix C]. See 
Appendix A.2 for a complete description of how all coefficients Ilan(?’T0 
and Zani'i~,r'), needed in the present work, are related to 0„_i(x)- 


6 An overview of the discretization 


Our Fourier-Nystrom discretization scheme is very similar to the scheme 
used in [20]. That scheme, in turn, builds on the schemes developed in m 
E2| in a pure Helmholtz setting. This section only gives a brief review. 

The FFT operations are, basically, controlled by two problem dependent 
integers N and Acorn with N > Neon, as follows: we use 2N + 1 equispaced 
points in the azimuthal Fourier transforms of kernels when r and r' are far 
from each other, we use 2Acon+l—u terms in the truncated convolutions m 
Equation (27)], and we use 2Acon + 1 equispaced points in the azimuthal 
Fourier transforms of smooth parts of kernels when r and r' are close. The 
value of N is chosen in an ad hoc manner. The value of Neon is determined 
by the decay of the Fourier coefficients of the test function 


9(0) 


sin(2A:pmaxSin(6>/2)) 

sin(0/2) 


(47) 
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where pmax is the largest value of p on 7 , so that Neon is the number of 
Fourier coefficients pn, n = 1,2,..., with Pnlgo > lOcmach- When r and 
r' are in the vicinity of corners, however, we may use smaller values of 
Noon- See, further. Section [7?^ We note, but do not generally exploit, that 
more elaborate adaptivity in the control of the FFT operations can lead to 
substantial computational savings. 

Our Nystrom discretization of (31), (32), (33) and (35) relies on an un¬ 
derlying panel-based 16-point Gauss-Legendre quadrature with a mesh of 
npnn quadrature panels on 7 . The IGupan discretization points play the 
role of both target points r* and source points Vj. The underlying quadra¬ 
ture is used in a conventional way when r* and Vj are far from each other. 
When Vi and rj are close and convolution is used, see Section an explicit 
kernel-split special quadrature is activated. Analytical information about 
the singularities in Kan{r,r') and San{r,r') is exploited in the construction 
of 16th order accurate weight corrections, computed on the fly. As to some 
extent compensate for the loss of convergence order that comes with the spe¬ 
cial quadrature, a procedure of temporary mesh refinement (upsampling) is 
adopted. See [18] for additional information on quadrature construction and 
upsampling. 

It is worth emphasizing that all iFo„(r, r') and San{r,r') contain some 
sort of singularities at r = r' and that these singularities are inherited by 
the corresponding Danir,r') and Zan{f,r') listed in Appendix A.2 The 
singularities are generally of logarithmic type, native to 0 ^_i (x), with the 
exceptions that the coefficients in Appendix A.2 that contain the function 


d{v) 


V ■ {r — r') 



v = T,iy,p,z, 


(48) 


may exhibit Cauchy-type singularities as r approaches r' and those coef¬ 
ficients that are proportional to yinix) have logarithmic-type singularities 
only in their first derivatives. The quadratures constructed in [18( I19| cover 
all these situations. 


7 Recursively compressed inverse preconditioning 

Spectral properties of integral operators in boundary integral equations are 
often sensitive to boundary smoothness. The very nature of solutions may 
be affected by a change in smoothness as may the performance of numeri¬ 
cal solvers. For example, the introduction of boundary singularities such as 
edges and corners can induce diverging asymptotics in layer densities. In¬ 
tense and costly mesh refinement is then needed for resolution, which may 
lead to the loss of stability. See [Hj for a review of recently developed 
numerical techniques to deal with this problem. 
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RCIP is one of the techniques discussed in [M]. It can be viewed as 
a general method to enhance the performance of panel-based Nystrom dis¬ 
cretization schemes. Roughly speaking, for Fredholm second kind integral 
equations, RCIP obtains solutions on piecewise smooth curves with the same 
accuracy and at about the same cost as with which solutions normally are 
obtained on smooth curves. The RCIP method originated in 2008 in the con¬ 
text of solving Laplace’s equation in piecewise smooth planar domains [21] 
and has since then been improved and extended as to apply to a variety of 
boundary value problems. 

A comprehensive description of the RCIP method is given in the tuto¬ 
rial m- This section first gives a brief summary and then focuses on some 
details particular to the MFIE system ( |31|) a nd (32) and to the converted 
normalization integral (35) with (36) and (37). 


7.1 Basics of the RCIP method 

Assume the following: we have an integral representation of a field U{r), 
r G ]R^\ 7 , in terms of an unknown layer density (7{r) on a piecewise smooth 
boundary 7 . On 7 there are a number ncm of corners with vertices 7 ^, j = 
1,... ,ncrn- The integral representation together with boundary conditions 
lead to a Fredholm second kind integral equation 


{I + K)a{r) = g{r), rGy, (49) 

where K is an integral operator with kernel K{r,r') on 7 . The operator K 
is compact away from the corners. The function g{r) is a right hand side 
with the same smoothness properties as 7 . We also assume that there is a 
relatively coarse mesh with Upan coarse quadrature panels of approximately 
equal length constructed on 7 . The purpose of the coarse mesh is to allow 
for a discretization that resolves g'(r), 7 , and K(r, r') for r far away from r'. 
We split the kernel 

K{r,r') = K*{r,r') + K°{r,r') (50) 


in such a way that K*{r, r') is zero except for when r and r' both lie within 
a distance of two coarse quadrature panels from the same 7 j. In this latter 
case K°{r,r') is zero. The kernel splitting (50) corresponds to an operator 
splitting 

K = K* + K\ (51) 

where K° is a compact operator. The variable substitution 

a{r) = {I + K*)~^ d{r) (52) 

lets us rewrite ( |4^ as a right preconditioned integral equation 

a{r) + K°{I + K*)~^a{r) = g{r), n E 7 , (53) 
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where the operator composition K°{I + K*) ^ is compact. 


The functions fj(r) and g{r) and the operator K° in (53) should be easy 
to discretize and to resolve on the coarse mesh. Only the inverse (/ + K*)~^ 
needs a hne mesh for its resolution. This fine mesh is constructed from the 
coarse mesh by, for each vertex 7 ^, choosing a number Usubj and letting the 
panels closest to jj be Usuhj times repeatedly subdivided. The size of Usuhj 
is determined by the behavior of (T{r) close to 7 ^ and by application-specific 


needs for resolution. The discretization of (53) can then be carried out as 


(Icoa T ^coa^) ^coa — Scoa ; 

where the compressed weighted inverse matrix R is given by 

R = P^ (I + K*)f:>. 


(54) 


(55) 


In (54) and (55) the subscript “coa” indicates a grid on the coarse mesh, the 


subscript “fin” indicates a grid on the fine mesh, the prolongation matrix P 
performs polynomial interpolation from the coarse grid to the fine grid and 
P^ is the transpose of a weighted prolongation matrix such that 


P^P — Icoa • 


(56) 


See [T71 Sections 4 and 5] for details. With 16-point composite quadrature 
the system size in ( |54[ ) is Ihripan x Ihripan- The matrix R differs from the 
identity matrix by having ricrn diagonal blocks R(-^\ j = 1,, ricm, of size 
64 X 64. 

Once (54) is solved for cfcoa, a discrete weight-corrected version of the 


original layer density is obtained from 


= R^r 


(57) 


The density <Tcoa, together with the composite quadrature, can be used for 
the accurate discretization of any integral on 7 involving a{r) and piecewise 
smooth functions. Furthermore, the field U{r) can now be recovered in 
those parts of the computational domain that lie away from the vertices 7 ^ 
using <Tcoa together with the quadratures of [ElDra] in a discretization of 
the integral representation for U{r). 

Note that in (54), the need for resolution in corners is not visible. The 
transformed layer density a-coa should be as easy to solve for as the original 


layer density in a discretization of (49) on a smooth 7 . All computational 


difficulties are gathered in the construction of the matrix R. 

There are 64 -|- 32nsubj discretization points on the fine grid within a 
distance of two coarse panels from the vertex 7 ^. Judging from an inspection 
of (55) it seems as if computing the matrix block should be an expensive 
and also unstable undertaking for large nsubj- Fortunately, R^-T can be 
computed via a fast and stable recursion which relies on hierarchies of small 
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(i) 

local nested grids around 7 ^ and produces hierarchies of matrices R- 7 i = 
l,...,ngubj, where the last matrix is equal to This fast recursion 

enables the computation of at a cost only proportional to Usubj- This 
is the power of the RCIP method. 

The fast recursion for R can also be run backwards, acting on d-coa,, 
for the purpose of reconstructing the solution crgn to a straight-forward 
discretization of (49) on the hne mesh 

(Ifin T Rfin) ^fin — Sfin • 


By this, one sees that the information contained in cTcoa; together with the 
Rp\ is the same as the information contained in crgn. A partial reconstruc¬ 
tion of cTfin is needed when U{r) is to be evaluated close to the vertices 77 
See [T71 Section 9] for a description of the reconstruction procedure. 


7.2 Details particular to the MFIE system 


The MFIE system (|31[) and (32), which on block operator form reads 


I - 2K^n 2ikS5n -2kSen 


Qsn{r) 


■ 0 ■ 

0 / -1- Kin i-R2n 



= 

0 

0 iRsn I + Kin 




0 


(59) 


has a more complicated appearance than the model equation (49). The 
operator corresponding to K* in (51), upon discretization, no longer yields 
a block diagonal matrix but a sparse block matrix where each vertex 7 ^ 
generates seven non-zero 64 x 64 blocks. In practice, this poses no problems 
for RCIP. Equation (54) still holds with 


Scoa — 0 ) 



Qsn 


r - 2 K° 

vn 

2iA:S5„ 

-2kSln ' 

^coa — 

Jrn 

K° = 

’ -^^coa 

0 


iKL 


^On 

coa 

0 

iKL 

1_ 


The compressed weighted inverse matrix 


of (55) is 


R = 


Pw 0 0 

0 P\Y 0 

0 0 Pw 


-1 

hH 

1 

to 

2ikS^n 

-2kSln ■ 

-1 

P 

0 

0 ■ 

0 




0 

p 

0 

0 

iKL 

I + K^n J 

fin 

0 

0 

p 


(60) 


(61) 


and has size 48npanX48repan- It can be permuted as to differ from the identity 
matrix by Ucm diagonal blocks R^l-* with 7 x 64 x 64 non-zero entries each. 
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When computing the via the fast recursion we take advantage of the 
sparsity structure in (61). We also allow for integers Wonj; controlling FFT 


operations and convolutions close to the vertices 7 ^, that may be smaller 
than the Won used on the coarse grid. 


These Wonj are determined as in 


Section ^ but with pmax of (47) replaced with the largest value of p on 7 


within a distance of two coarse panels from 7 j. 


7.3 Resolving the normalization integral 

The RCIP method provides a tool for the fast and accurate solution of 


the MFIE eigensystem (59) within the framework of our Fourier-Nystrom 
scheme. The discretized equation (54) with (60) and (61) can be used both 


to find eigenwavenumbers and to hnd the corresponding discrete transformed 

eigenvectors, that is, non-trivial solutions. The eigenvectors can, in turn and 

f?) 

together with the matrices R^ 7 be used to reconstruct the discrete densities 
Qsm Jrnj uud jQn on the fine grid. This is enough to allow for the accurate 
evaluation of non-normalized electric eigenfields in the entire computational 
domain, but not enough to allow for the accurate evaluation of normalized 
eigenfields. 


The converted normalization integral (35) with (36) and (37) requires 


that the Fourier coefficients of 'l'(r) and A(r), and their derivatives with 
respect to r and , are sufficiently resolved on the fine grid so that their 
various inner products and squared moduli can be accurately integrated 
along 7 . For a prescribed overall accuracy and for densities Qsnir) or Jenir) 
that diverge at corner vertices, this poses tougher requirements on panel 
refinement than merely demanding that the densities are sufficiently resolved 
as to be accurately integrated against piecewise smooth functions. Mappings 


from partially reconstructed values of 


and Jgri, on the fine grid to 


values of these sought Fourier coefficients on the fine grid can be performed 
via hierarchical matrix-vector multiplications similar to, but simpler than, 
those used for the reconstruction of Jrn, and Jg^ themselves. This 
procedure uses hierarchies of small matrices corresponding to the evaluation 


of all modal operators present in (40), (41) and (42) on the small local nested 


grids mentioned in Section 7.1 


8 Physical fields and edge singularities 


This section relates complex valued electric eigenfields in of the form (34) 


to the physical time-domain standing wave fields that are excited and mea¬ 
sured in real life experiments. To facilitate the interpretation of, so called, 
field maps we also review the leading order asymptotic behavior of electric 
fields and surface charge and current densities close to edges. 
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8.1 Physical time-domain fields. 

Every eigenvalue k'^ of the system §, § , and ([^, not belonging to the 
mode n = 0, is degenerate and typically corresponds to a two-dimensional 
subspace of electric eigenfields. Such an eigenspace can be spanned by two 
orthonormal eigenhelds En{r) and of the form (34), constructed 

via (33) from the Fourier coefficients gsn{r), Jrnir), £>s(-n)(^)) 

JT(-n)\'f')i that are non-trivial solutions to (31) and (32) at eigen- 

wavenumber k and normalized with (|^. 

The normalized Fourier coefficients are, in turn, unique only up to a 
constant factor of modulus one. In our implementation we choose these 
factors such that Jmir) is real and even in n, that is, = Jrni'f')- 

Equations (31), (p^, (p3|) and the formulas in Appendix A.2|then imply the 


following: Jeri{f) is imaginary and odd in n; Epn{r), E^nir), and Qsni'r) are 
imaginary and even in n; and Egnir) is real and odd in n. 

Complex valued standing waves are formed by linear combinations of 
Fourier coefficients as 


4^r!{r) = 


Jr(-n){r)e = anJrnir) cos n6 , 
4n{r) = -iY(«^rn(?')e‘"® “ <^T(-n)= ttn Jrnir) smn6 , 
where 


(62) 


ao = l/V^, a„ = l/V7r, n = l,2 ,..., 

and superscript (e) and (o) denote “even” and “odd”. The prefactor —i in 
the second equation of (62) is to make both surface currents real. 

The physical time-domain currents in the r-direction are now obtained 
from 


(r, t) = Re{ (r)e = anJmir) cos n9 cos uit , 

= Re{ = anJT-n(?') sin n0 cos wt. 


(63) 


Expressions for the other physical time-domain quantities are uniquely com¬ 
posed in the same manner as 

n(e)/ 


= -ianQsn{r) COS n9 sm Ujt, 
= ianJenir) sm n9 cos ut, 
E^^{r,t) = —ia„Fp„(r) cosn^sincjf , 
= anEgn{r) S\nn9 sinut, 
E4}{r,t) = —ianEzn{r) cos n9 sin cot, 


(64) 
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and 


= —ianQsnir) sinn9smujt, 
—ianJeni^) cos nO cosut, 
—ianEpn{r) sin nO sincji , 
—anEon (r) cos n9 sin cot , 
E^^{r,t) = —ianEznir) smn9 sincot. 




E. 


pn 


(65) 


All components of the physical electric eigenheld and its surface charge den¬ 
sity are in phase with respect to time, but 90 degrees out of phase with the 
surface current density. In the numerical examples of Section we present 
field maps in the xz-plane {9 = 0, tt) of the imaginary part of Epn{r)e^^^ and 
Ezn{r)E'^^, and the real part of E 0 n{r)e^^^. 


8.2 Asymptotic behavior at edges 

Let a corner with vertex 7 ^ have an (inner) opening angle of aj. Let be 
the tangential distance from a point r G 7 to jj and let ^ be the Euclidean 
distance from a point r G AuL to 7 ^. For n > 0, we have the general leading 
behaviors close to 7 ^ 


Qsn{r) ~ ^ 

Epn{r) ~ , 


Een{r) ^ 1 , 


Jen{r)-^e^ ' 

Ezn{r) ~ 


where 


Pj = 


TT 


27r — an 


( 66 ) 


(67) 


See [S] and references therein. The asymptotics of the Fourier coefficients 
of T(r) and A(r), and their derivatives with respect to r and are more 
complicated. 


For the special case of n = 0 one can show that ( 66 ) holds with some 


sharpening: the even fields of (64) have Jeoir) = 0 and EgQ{r) = 0; the odd 


fields of (65) have £»so(^) = 0, Eoii") = 0, and Epo{r) = E^oir) = 0. 


(32), and (33) is implemented in 


9 Numerical examples 

Our Fourier-Nystrom scheme for 
Matlab and executed on a workstation equipped with an Intel Core i7- 
3930K CPU and 64 GB of memory. The weight corrected densities 
Jrn, and Jdn on the coarse grid are obtained from (54) with (60) and (61). 
The densities Jrn, and Jgn on the fine grid are obtained with recon¬ 
struction m Section 9]. To enforce ([^ we normalize the densities with the 
value of 11 Fin 11 obtained from their insertion in a discretized version of (35). 
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The Matlab implementation is standard and relies on built-in functions. 
No particular attempts are made at optimizing the code for speed, except 
for the use of a few parfor-loops (which execute in parallel). Great care 
has gone into obtaining intermediate quantities to high accuracy and to 
resolve modal integral operators in corners. The solution time quoted in the 
examples below refers to wall-clock time from when an eigenwavenumber is 
known and until the normalized densities J™, and Jon are obtained on 
the fine grid. 


9.1 Search for eigenwavenumbers 


Eigenwavenumbers are determined using a separate, slimmed down, code 
that is cleared from matrices and panel refinement only needed for the nor¬ 
malization. In what follows we let knj be the jth smallest eigenwavenumber 
for azimuthal index n. Our search algorithm for kn,j^ with n fixed, is a mod¬ 
ification of the “standard published method” described in |3l Appendix B]. 
The standard method is to search along the fc-axis for (near) zeros of the 
lowest singular value s{k) of an appropriate discretized system matrix B(fe). 
Successive parabolic interpolation, which has convergence rate Tc ~ 1.324, 
is applied to s‘^{k) and is safeguarded by the empirical observation that the 
slope of s{k) appears to have a domain dependent upper bound Cg of size 
0(1). In our setting, B(A:) can be taken as the lower right 32n 
part of -|- K°q 3 ^ with R and from (61) and (60). 


pan ^ 32?T.pan 


We modify the standard method, described above, by searching for zeros 
of the smallest eigenvalue X{k) of B(A;) rather than for (near) zeros of the 
smallest singular value. We then replace successive parabolic interpolation 
applied to s‘^{k) with the secant method applied to |A(A:)| sgn(arg(A(A:))). 
The slope of the function |A(A;)| sgn(arg(A(/c))) also appears to have a domain 
dependent bound of size 0 ( 1 ), which we denote Cx and use for safeguarding. 
The secant method has convergence rate rc ~ 1.618. When Cx is chosen 
correctly and for a fixed n, our modified search algorithm finds all knj in 
a given fe-interval typically needing between four and eight iterations per 
eigenwavenumber found. 


9.2 Comparison with solution in semi-analytic form 

The codes are first verified for V being the intersection of a cone with half 
opening angle of one radian and a spherical shell with outer radius one and 
inner radius 0.5. The generating curve 7 is parameterized as 

( 0.5(sin(2t-F2),cos(2t-h2)) , tG[-l,-0.5], 

^(^) = S (t + 1) (sin(l),cos(l)) , tG [-0.5,0], ( 68 ) 

[ (sin(l — t),cos(l — f)) , t G [ 0 , 1 ] , 

and has non-reentrant corners with vertices at 71 = 0.5 (sin(l), cos(l)) and 
72 = (sin(l),cos(l)). This cavity is an excellent test geometry since, while 
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Eigenfield E e'” at /(=13.724219259476561 
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Figure 2: The electric eigenfield at /ci_io = 13.724219259476561 for the cav¬ 
ity with 7 as in (68). Left: field maps in the xz-plane {6 = 0, vr) of (a) 
Im {£pi(r)e'®}, (c) Re {£' 0 i(r)e‘®} , and (e) Im {£ 2 i(r)e^®}. The null-field 
outside the cavity is omitted. Right: (b), (d), and (f) show log^g of the point- 
wise absolute difference between the field maps of our scheme and those from 
the semi-analytic solution. See Section |^for an explanation of how field maps 
relate to physical time-domain fields. 


not trivial, it allows for semi-analytic solutions to the original system (§,(©, 
and (©. The normalized eigenfields En{r) are expressed in regular and 
irregular spherical vector waves, see [HI Section 2.2], that are modified to 
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satisfy the boundary condition Q. The major modifications are that the 
associated Legendre functions P”( 2 ;/|r|) and spherical Bessel and Neumann 
functions, ju{k\r\) and yu{k\r\) in the vector waves, have indices v and 
wavenumbers k that are solutions to transcendental equations. 

Figure [ 1 (a), He), and He) show field maps of Eai{r)e^^^ a = p,0,z, at 
eigenwavenumber = 13.724219259476561 constructed from Fourier co¬ 
efficients Eai{r) produced by our codes. The eigenwavenumber corresponds 
to about 3.7 wavelengths across the generalized diameter of V. The Fourier 
coefficients are obtained with Upan = 28 quadrature panels, corresponding 
to 448 discretization points on 7 , and they are evaluated at 245000 points 
on a Cartesian grid in the rectangle 


Oi = {r G : 0 < p < 0.95, -0.45 < z < 1.45} . 


(69) 


Only coefficients with r ^ A are actually used in Figure [^a) , [^c) , and[^e). 
The FFT operations are controlled by the integers N = 141, iVconi = 24, and 
.^con 2 = -^con = 30, See Sectionsand 7.2 The densities Qso{r), Jroir), and 
Jeo{r) are bounded and the number of panel subdivisions used by the RCIP 
method for resolution close to 71 and 72 is chosen as ngubi = f^sub 2 = 30, see 


Section 7.1 The solution time is around 13 seconds and the time required to 
evaluate the coefficient vector {Epi{r),Egi{r),Ezi{r)) is, on average, 0.002 
seconds per point r. 

Figure H^)’ H*^)’ 0^) show log^o of absolute difference be¬ 

tween the field maps produced by our codes and the field maps obtained 
from the semi-analytic solution. Here all coefficients with r G fli are used. 
When obtaining the semi-analytic solution, the eigenwavenumber is evalu¬ 
ated to machine precision and the eigenfields to almost machine precision 
using a combination of Matlab with extended precision and Maple. The 
semi-analytic solution at points r outside A U 7 is taken as Eai{r) = 0, 
compare ( p^ . One can conclude that, in this example, our codes give co¬ 
efficients Eai{r) that are pointwise accurate to at least 13 digits and an 
eigenwavenumber that is accurate to machine precision. 


9.3 The one cell elliptic cavity 

Our remaining numerical examples pertain to the cavity depicted in Figurej^ 
which, in particle accelerator terminology, is known as a one cell elliptic 
cavity. The generating curve 7 is parameterized as 


r{t) 


(vr-h t,-1 - 7r/4) , 
( 7 r/ 4 ,-l-k 7 r/ 2 -Ft) , 

< (vr/d-|-cos(t), sin(t)) , 
(vr/d, 1 - 7 r/ 2 -ht) , 

^ (vr-t, l-hvr/d) , 


t G [—vr, — 37 r/ 4 ] , 
t G [— 37 r/ 4 , —vr/2] , 
t G [- 7 r/ 2 , 7 r/ 2 ] , 
t G i7r/2,37r/4], 
t G [ 37 r/ 4 , vr] , 


(70) 
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and has corner vertices at 71 = (vr/d, —1 — vr/4), 72 = (vr/d, —1), 73 = 
(vr/d, 1), and 74 = (vr/d, 1 + vr/d). The corners at 72 and 73 are reen¬ 
trant. The number of panel subdivisions in the RCIP method is chosen as 
^T-subi = n-suhA = 30 and nsub 2 = ^T-subs = 140 in all examples. The integer N, 
controlling the FFT operations when r and r' are far from each other, is by 
default chosen as N = dnpan + n. The estimated pointwise absolute error in 
a given computed field map is based on a comparison with a more resolved 
map obtained with 50 per cent more quadrature panels on 7 . Fourier co¬ 
efficients are evaluated on Cartesian grids in rectangles, most often chosen 
as 

= {r G : 0 < p < 2,-2 < z < 2 } . (71) 

Superconducting elliptic cavities are common in linear accelerators for 
protons. They are to be used in a projected Superconducting Proton Linear 
accelerator (SPL) at CERN and in the linear accelerator for the European 
Spallation Source (ESS) that is currently under construction in Lund, Swe¬ 
den. In the design of elliptic cavities it is important to determine several 
quantities to high accuracy. For the fundamental eigenfield that accelerates 
the protons, one needs to evaluate the maximum normalized electric and 
magnetic fields on the surface and the normalized electric field on the sym¬ 
metry axis. The eigenfields with azimuthal indices n = 0 and n = 1 have 
non-zero held components on the symmetry axis which cause them to in¬ 
teract with the beam of particles. A large number of these eigenhelds have 
to be determined in order to assess their effect on the beam. The three 
numerical examples we present for the elliptic cavity have n = 0 and n = 1 
and eigenwavenumbers that are relevant for particle accelerators. 

9.4 The fundamental mode 

The fundamental electric eigenheld is the eigenheld with the lowest eigen- 
wavenumber. For the elliptic cavity with 7 as in ( |70| ) it has eigenwavenum- 
ber ko^i = 1.5631689906935362, corresponding to 0.97 wavelengths across 
the generalized diameter of V. Figure j^a) and i c) show field maps of 
£'po(^) and Ezoi^r) as computed with our scheme. The map of Eeo(r) is zero 
and therefore omitted. The Fourier coefficients are obtained with npan = 32 
quadrature panels, corresponding to 512 discretization points on 7 , and they 
are evaluated at 245000 points on a grid in fin- The FFT operations, for r 
and r' close, are controlled by A'coni = A"con 4 = 12 , A^con 2 = A^conS = 14 and 
Aeon = 16. The solution time is around 16 seconds and the time required to 
evaluate the coefficient vector {Epo{r), Ezo{r)) is, on average, 0.003 seconds 
per point r. 

Figure [^b) and [^d) show logj^Q of the estimated pointwise absolute 
error in Figure l^a) and[^ c). At this low eigenwavenumber the estimated 
accuracy is quite exceptional. The solver delivers 15 accurate digits except 
at points close to 7 . 


22 



Figure 3: The fundamental electric eigenfield for the elliptic cavity with 7 
as in (70). Left: field maps in the xz-plane of (a) Im{F^po(^)} and (c) 
Im{F^^o(^)}- Right: log^Q of estimated pointwise absolute error. 


Note that Ezo{r) is strong along the symmetry axis. This explains why 
the fundamental mode is used for acceleration of charged particles. At the 
vertices of the reentrant corners both Epo{r) and Ezo{r) diverge as 
see Section |8.2[ In the design of de facto elliptic cavities in accelerators, 


sharp reentrant cell- and iris edges are avoided. On the other hand, there 
are sharp reentrant edges where the beam pipe is attached to the cavity and 
it is therefore important that a solver can handle all sorts of sharp edges. 
We have omitted the beam pipe in order to keep the model simple. 


9.5 A convergence study 

The next example is the eigenfield with ^ 1^9923 = 120.2309391499240. De¬ 
spite a generalized diameter of the elliptic cavity that now corresponds 
to around 75 wavelengths, our solver maintains 16th order convergence 
and high achievable accuracy, as seen in the left image of Figure The 
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estimated average absolute error in A 


Eigenfields of E at /(=120.2309391499240 Behavior of layer densities close to y 



Figure 4: Left: convergence of £'pi(r), E 0 i{r), and Ezi{r) with 7 given by® 
and at eigenwavenumber A:i ,9928 = 120.2309391499240. The estimated aver¬ 
age pointwise absolute error in A has converged to less than 10“^^ at 1920 
discretization points on 7 , corresponding to 16 points per wavelength along 7 . 
Right: behavior of fei(r), Jriir), and J 6 ii(r) close to 72 . 


FFT operations are in this study controlled by = max{300,4npan + n}, 
IVconi = lVcon4 = 127, Alcon2 = NconS = 137, and iVcon = 259. The Fourier 
coefficients are evaluated at 45000 points on a grid in fin and have con¬ 
verged to more than 12 digits already at 16 points per wavelength along 7 , 
which is marginally better than in a similar study for the eigenfield with 
^ 1,2460 = 60.21392380136615 (not shown). We conclude that there are no 
signs of any pollution effect [ 2 ] at these wavenumbers. 

The right image of Figure [^reveals that at the reentrant corners both 
Qsi{r) and J 0 i{r) diverge as , whereas Jriir) is bounded. These asymp¬ 


totics are in accordance with ( 66 ). 

Figure shows fully converged field maps of Eai{r)e^^, a = p,6,z, ob¬ 
tained with 2816 discretization points on 7 , along with log^Q of estimated 
pointwise absolute error. The solution time is around 240 seconds and the 
time required to evaluate the coefficient vector {Epi{r), E0i{r), E^iir)) is, 
on average, 0.04 seconds per point r. One can see, in the left images, that 
the eigenfield is concentrated to a region close to the symmetry axis. This 
is typical for eigenfields with small n and large k. 

It follows from (10) that the normal component of the coefficient vec¬ 
tor has the same (singular) behavior as Qsiir) along 7 . The amplitude of 
a singularity in ^)si(^’) in a reentrant corner is often small at large eigen- 
wavenumbers. As seen in the right image of Figure]^ it may become visible 
first at a distance from a corner vertex that is less than one thousandth of 
the total arclength. Although the images of Figure use 245000 evaluation 
points on the grid in Qu, this is not enough to detect the singularities in the 
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Figure 5: The electric eigenfield for the elliptic cavity with 7 as in (70) and 
with A:i ^9928 = 120.2309391499240. Left: field maps in the xz-plane of (a) 
Im {£'pi(r)e'®}, (c) Re {£' 0 i(r)e‘®}, and (e) Im {£i 2 i(r)e'®}. Right: log^o 
estimated pointwise absolute error. 


field maps. This underscores the importance of being able to zoom regions 
where singularities might appear in order to determine their amplitudes. 
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Eigenfield at k=3^ .65910852052012 


Log^Q of estimated absolute error in 



Figure 6: The electric eigenfield for the elliptic cavity with 7 as in (70) and 
with A:o ,662 = 31.65910852052012. Left: (a) the field map of Im{£'j;o(?’)} in 
the xz-plane, (c) is a ten times and (e) a 100 times magnification in the vicinity 
of the corner vertex 73. Right: log^g of estimated pointwise absolute error. 
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9.6 Corner zoom 


Our last example is the eigenfield with A:o 662 = 31.65910852052012. This 


eigenheld is even according to the classification of Section 8.1 and has 
E 0 o{r) = 0 in agreement with Section [8^ Figure [^a) shows a held map of 
Ezo{r) obtained with 768 discretization points on 7, iVconi = .^con4 = 48, 
I^con 2 = = 57, and Neon = 85 and evaluated at 245000 points on a 

grid in fin- The solution time is around 37 seconds. 

Figure l^c) and i e) explore the held map of E 2 Q{r) when the region 
around the corner vertex 73 is magnihed hrst ten times and then 100 times. 
There are still 768 discretization points on 7, but held evaluations now take 
place at 490000 points on grids in the squares 


= <^ r e 


TT ^ TT 

— — 5 < p < - 
4 “ 4 


+ (5,1 — (5< 2:<1 + (5| 


(72) 


where 6 = 0.2 or (5 = 0.02. At 73 (and at 72) the held map of E^oir) exhibits 
the singularity given by (66). The right images of Figure]^ show log^Q of the 
estimated pointwise absolute error. 

This example emphasizes that high accuracy is vital for the detection 
of singular helds. A less accurate solver might neglect strong local electric 
helds that can lead to serious electric discharges in a cavity. 


10 Conclusion and outlook 

We have presented a competitive solver for the determination of normalized 
electric eigenhelds in axially symmetric microwave cavities with piecewise 
smooth PEC surfaces. The solver is based on the following key elements: 
the interior magnetic held integral equation, a charge integral equation, a 
surface integral for normalization, a high-order convergent Fourier-Nystrom 
discretization scheme, on-the-hy computation of singular and nearly singular 
quadrature rules, and access to high-order surface information. 

While our solver determines eigenhelds with extraordinary accuracy for 
a large range of eigenwavenumbers, we observe that the development of effi¬ 
cient high-order Nystrom schemes for time-harmonic boundary value prob¬ 
lems in piecewise smooth three-dimensional domains is an active research 
held. See [7] for a recent example of a fast solver for the integral equations 
which model low-frequency acoustic scattering from curved surfaces. 

The computational needs in accelerator technology are extensive and our 
solver must be equipped with additional features to become a truly versa¬ 
tile tool. Our next step is to allow for sources, modeling a pulsed beam 
of particles, inside cavities. This makes it possible to evaluate wakehelds 
generated by beams in accelerators. We foresee that our solver can be used 
for benchmarking. It will be able to evaluate high-frequency parts of wake- 
held spectra that other solvers cannot reach and by that it becomes an 
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important complement to state-of-the-art software such as the CSX Parti¬ 
cle Studio Wakefield Solver. We also anticipate other improvements. The 
present Matlab implementation places emphasis on high achievable accu¬ 
racy and on rapid convergence with respect to the degrees of freedom. For 
this, modal integral operators are upsampled in a somewhat crude and costly 
way. Better adaptivity in this process and in the FFT operations will lead 
to increased execution speed. 

Nano-optics is another application that involves high-frequency electro¬ 
magnetic eigenfields. Here structures that support whispering gallery modes 
(WGMs) are of great interest [27]. The WGMs have large eigenwavenumbers 
and large azimuthal indices and the numerical examples presented in [20] 
indicate that our solver is ideal for their evaluation. When the WGM struc¬ 
tures are dielectric bodies of revolution, the solver needs to be adapted to 
a similar set of integral equations as in [^. The WGM structures used in 
nano-optics are designed to have very high Quality factors (Q-factors), and 
the bandwidth of such structures is very small since the Q-factor equals the 
frequency-to-bandwidth ratio. There are reported Q-factors as high as 10^^ 
m with a relative bandwidth of 10 In principle, a solver that delivers 
eleven digit accuracy is needed to design such a structure if a WGM of the 
structure is to be excited by a laser with a given wavelength. 
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A Explicit expressions for kernels 

The various double- and single-layer type operators and Sa used in 
this paper are defined by their corresponding static kernels Da(r,r') and 
Za(r,r') via (18) with (19) and (20) with (21). This appendix collects 
explicit expressions for all static kernels along with analytic expressions for 
their Fourier coefficients Dan{'r,r') and Zan{f,r'). The numbering of the 
operators and kernels is compatible with the numbering used in m- Note 
that for azimuthal index n = 0 the modal operators K 2 n, K^n, Ri 2 n, Ki^n, 
K\Qni K 2 Qni I^ 2 ini S 2 ni 5*3^1, Sgn, and 5*1372 are zero. 


A.l Static kernels 

The static kernels are expressed in terms of the azimuthal angle 9 and vectors 
and points in the plane defined by 0 = 0, see Section [2.1] The abbreviations 

1 / ■ {r — r') = v ■ {r — r') + Upp'{l — cos(9 — 6 ')) , 
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T ■ {r — r') = T ■ {r — r') + Vzp'i^ — cos(0 — 9 ')), 
f — r'l = {p^ + p'^ — 2pp' cos{6 — o') + {z — z')"^) ^ , 


^oir,r') = 


1 


mr,r') = 


1 


47r|r —r'l’ ’ 47r|r —r'P’ 

are used. The static double-layer type kernels are 


Di,{r, r') = -u {r - r')^{r, r'), 

Di{r, r') = 2 {u'^p - {u' ■ r' - v'^z) cos{e - 9')) r '), 

D 2 {r, r') = — 21 ( 2 ; — z') sin(0 — 0')$Q(r, r '), 

T) 3 (r, r') = 2\{v'^v • r — • r') sin(6* — 6*')$o(r, r '), 

Di{r, r') = -2 {vpp - {v ■ r - v^z ) cos(6» - 9')) $o(r, r '), 
Dii(r, r') = [p- p cos{9 - 9')) io(r, r'), 

Duir, r') = -ip sin(6' - 6'')d>o(r-, r'), 

Di 3 {r,r') = {z - z')^{r,r'), 

Di 4 {r, r') = —T • (r — r')v'^ cos{9 — 6*')4>o(r, r '), 
i4i5(r, r') = It • (r — r') sin(0 — 0')$o(r, r'), 

DiQ{r, r') = —IT • (r — r')v'^ sin(0 — 0')$Q(r, r '), 

Dij{r, r') = —T • (r — r') cos{9 — 0')io(r, r '), 

Dis{r, r') = T -{r - r')v'p^l{r, r '), 

Di^{r, r') = —u • (r — r')v'^ cos(0 — 0')$o(r, r'), 

D 2 o{r, r') = iv • {r — r') sin(0 — 0')$Q(r, r '), 

D 2 i{r, r') = —\u ■ (r — r')u'^ sin(0 — 6*')'lQ(r, r '), 

D 22 {r, r') = —u • (r — r') cos(9 — 0')4>o(r, r'), 

J^ 23 (r, r') = ly- (r- r')iy'p^^(r, r '). 


The static single-layer type kernels are 


Z^{r,r') = 4'o(r,r') 

Zi{r,r') = {uziy'zCos{9 - 9') + Vpv'p) 4'o(r*,r'), 
Z 2 {r, r') = -ivz sin(0 - 6 »')^o(t‘, r '), 

Z^{r, r') = iu'^ sin(6' - 6'')$o(j', r '), 

Z 4 {r, r') = cos{9 — 9')^o{r, r '), 

Zbir, r') = [upu'^ cos(6' - O') - v^^'p) r'), 

ZQ{r, r') = -ivp sin(6» - 9')^Q{r, r'), 

Zj{r, r') = z/' cos(0 — 9')^o{r, r '), 

Zs(r,r') = isin(6' - 9')^o{r,r'), 
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Z9(r,r') = -Z3{r,r') , 
Zio{r,r') = Z4{r,r'), 
Zn{r,r') = -Up<^>oir,r '), 
Zi2{r,r') = Zj{r,r '), 
Zn{r,r') = -Z8{r,r') , 
Zi4(r,r') = Zn{r,r'). 


A.2 Fourier coefficients 


Our derivation of the Fourier coefficients of the static kernels is made in a 
similar manner as in Young, Hao, and Martinsson [321 Section 5.3]. The idea 
to expand the Green’s function for the Laplacian in the functions 0„_i(x) 
comes from Cohl and Tohline [lO]. We use the notation of Sections 2.1 and[^ 
with X as in (44) and 


r] = ■ 

V ■ {r — r') 


d{v) = 


\r — r 


./12 ’ 


'fin(x) = \ (53„(x) +n„_i 
^nix) = (x) - , 

£4x) = 51„(x) + 2(x-l)Vn(x), 
^n{x) = 2n(x - 1)£1„_ 1 (x) ■ 


The Fourier coefficients of the static double-layer type kernels are 


D^n{r, r') 
T>in(r-, r') 
D 2 n{r, r') 
Dsnir, r) 
D4n{r, r') 
Diin{r, r') 


V 
— 2rj 


d^^nix) - —^n(x) 
p 




pp 


(z — z') ^ 

PP ^ 

(iz'v • r — • r') ^ 


PP 


—2r] 

-V 


dimnix) + 


PP 


1 


d{p)^n{x) - -*fin(x) 
P 


30 



















Di2n{r,r) = --nn„_i(x), 

P 2 

Di3n{r,r') = -pd{z)d\n{x ), 

DlAn{r, r') =r]iyi d{T)2n{x) - ^ {^n{x) + ^n{x)) , 

Dl5nir, r') = P d{T)ttnix) - ^ (<^n(x) + ®n(x)) , 

L 2/9 J 

DlSnir, r') = -p d{T)^n{x) - ^ {^n{x) + Sn(x)) , 
Dnn{r, r') = p (i(r)£„(x) - ^ (£n(x) + 5Jt„(x)) , 
Di3n{r,r') =-pv'p d(r)9^„(x) - ^^n(x) , 

^19n(r', r') = 7? 12 ^ d{iy)£n{x) - ^ ('2n(x) + 5?l„(x)) , 
D 20 n{r, r) = p d{v)(tn{x) - ^ {^n{x) + Sn(x)) , 
D 21 n{r, r) = -p d{v)(tn{x) - ^ {^n{x) + Sn(x)) , 
D 22 n{r, r') = P d{u)£,n{x) - ^ (-C„(x) + 9?tn(x)) , 

D23n{r,r') = -pu d{iz)na{x) - —^nix) ■ 

I Pi 

The Fourier coefficients of the static single-layer type kernels are 

Z,n{r,r) = pQ^_i{x), 

Zin{r,r') = p i^piy'p£l^_i{x) , 

Z2n{r,r') = pu^'Snix ), 

^3n(7’,r') = -pv'^Ttnix) , 

ZAn{r,r’) = pTlnix), 

Z5n{r,r) = P Vpv'^Tlnix) - Vz^'pQn-\bd , 
^6n(7’,r') = pVpDnix) , 

^7n(?',r') = pv'^Tlnix) , 

Zsn{r,r') = -r?D„(x), 

^9n(r,r') = -Z3n{r,r') , 

^iOn(?'T0 = ZAn{r,r') , 

^iin(r-,r') = -py'p£i^_i{x ), 

^12n(?'T0 = ^ 7 n(r,r'), 
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Zi3n{r,r') = -Zsnir,r') 
Zun{r,r') = Ziin{r,r'). 
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